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ABSTRACT 


A steady, axially symmetric, incompressible, viscous flow past a rigid sphere is 
numerically simulated by using a numerical scheme, based on spectral methods. The 
equations have been reduced to two sets of nonlinear second order partial differential 
equations in terms of vorticity and stream function. The calculations have been carried out 
for Reynolds numbers, based on the sphere diameter, in the range 0.1 to 104. 

The numerical results have verified that there is excellent agreement with Stokes 
theory at very low Reynolds numbers. At moderate to intermediate Reynolds numbers 
there is good general agreement with available experimental data and flow visualization 
pictures. The Reynolds number at which separation occurs is estimated as 20. The 
approach to boundary-layer behavior with increasing Reynolds numbers is also verified by 


comparison with potential flow theory and analytical boundary-layer solution. 
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I. INTRODUCTION 


Numerical simulations in fluid dynamics have historically attracted considerable 
attention, through the use of various numerical schemes such as finite-difference and 
finite-element methods. In this study, the direct numerical simulation of a time-dependent, 
axially symmetric, incompressible viscous flow past a rigid sphere was studied by using a 
different numerical scheme based on spectral methods. Spectral methods have become 
increasingly popular in recent years, especially with the development of fast transform 
methods. It has been successfully applied to numerical weather predictions, numerical 
simulation of laminar/turbulent flows and other fundamental problems especially where 
high accuracy is desired for relatively lower memory requirements. 

The drag dependence on Reynolds number was studied and the results were 
compared with the previous results. The stream function wy and the vorticity w were 
visualized through the use of the numerical results and compared with the flow 
visualization pictures for various Reynolds numbers. The boundary-layer thickness on the 
surface of the sphere was also investigated through the use of the results of the 
simulation. The boundary-layer thickness on the surface of the sphere is calculated by 
introducing the power series expansions of the stream function into the boundary-layer 
equations. The results are valid for Re > © analytically. The minimum Reynolds number 
for which the boundary-layer assumption is valid, is investigated by comparing the 
boundary-layer thickness which is obtained from the numerical solution, with analytical 
results in the boundary-layer limit 

It is clear from the foregoing that the objective of this study is to carry out the 
numerical scheme based on the spectral methods, through the use of the vorticity-stream 
function form of the Navier-Stokes equations, on a time-dependent, axisymmetric, viscous 
flow about a rigid sphere. A Matlab program was developed to implement the foregoing 
numerical scheme. The results were compared to the previous studies which were based 


on different methods. However the details of the method itself were not compared to the 


other methods which were used in previous studies. The expectations are that the results 
will match with the previous results, improve our current ability to compute/predict the 
evolution of flow for a particular geometry and conditions, and, hopefully shed some light 


on the physics of flows which have been heretofore uncalculated. 


Il. BACKGROUND STUDIES 


The collocation method appears to have been first used by Slater (1934) and 
Kantorovic (1934) in specific applications. It was developed as a general method for 
solving ordinary differential equations. In 1938 Lanczos’ work showed for the first time 
that a proper choice of trial functions and distribution of collocation points is crucial to the 
accuracy of solution. This method was revived by Clenshaw (1957), Norton (1963) and 
Wright (1964). The application of Chebyshev polynomial expansions to the initial value 
problems are involved in these studies. 

The spectral collocation method was applied to partial differential equations for 
spatially periodic problems by Kreiss and Oliger (1972) (who called it the Fourier method) 
and Orszag (1972) (who termed it pseudospectral). These were the earliest applications of 
spectral collocation or pseudospectral method. This approach was very attractive because 
of its application to variable-coefficient and even non-linear problems. 

The Galerkin approach which depends on the same trial and test functions, was 
applied to a meteorological model by Silberman (1954). This was the first serious 
application of spectral methods to PDE’s. After Orszag (1969,1970) and Eliasen, 
Machenhauer and Rasmussen (1970) developed transform methods for evaluating 
convolution sums arising from quadratic non-linearity, spectral Galerkin methods only 
became practical for high resolution calculations of such non-linear problems. 

The first unifying mathematical assessment of the theory of spectral methods was 
covered in the monograph by Gottlieb and Orszag (1977). Since then, the theory has been 
extended to cover many different problems, such as variable-coefficient and non-linear 
problems. Applications in fluid dynamics were reviewed in the symposium proceedings 
edited by Voigt, Gottlieb and Hussaini (1984). Other introductory or review articles have 
appeared recently such as Mercier (1981), Hussaini, Salas and Zang (1985), Zang and 
Hussaini (1985), Deville (1984), Gottlieb (1985) and Hussaini and Zang (1987). 


It is well known that the background of this particular problem of flow over a 
sphere has a long history due to the fundamental interest in external flow over bluff 
bodies. The first theoretical results for steady, uniform flow past a rigid sphere were given 
by Proudman and Pearson (1957) and are valid for Re < 1. For higher Reynolds numbers 
computations have been conducted by Rimon and Cheng (1969), Le Clair, Hamielec and 
Pruppacher (1970) and more recently by Fornberg (1988). All these computations are 
based on finite-difference schemes. Dennis and Walker (1971), Oliver and Chung (1985) 
and Brabston and Keller (1975) used a different computational approach based on an 
expansion of Legendre functions and spherical harmonics. The results of these simulations 
are consistent with the experimental data for the drag forces, and have been used to 
develop improved formulae. Marcus and Tuckerman (1987) developed another technique 
based on the spectral methods and successfully applied it to the simulation of a flow 
between concentric rotating spheres. And more recently Chang and Maxey (1994) have 


used the same method to compute the oscillatory flow about a sphere. 


It. GOVERNING EQUATIONS 


A. DERIVATION OF THE EQUATIONS 


The well known Navier-Stokes equations of motion of a viscous, incompressible 


fluid are; 


du 1 ‘ 
>, go u (aly) 


V.u=0 (3.2) 


Define dimensionless variables 








nr . u . Vp . i 
i Oe —= V=- V 
6 a U pU° a 
where a: radius of the sphere 
U: free stream velocity 
p : density of the fluid 
v : kinematic viscosity of the fluid 
Introduce the dimensionless variables in Eq. (3.1) 
Uvdu U 1 U’ U 
——————+—_(1 Ve =—— ee ae 
a’ ot 7s i Op a ele ie o>) 
Divide Eq. (3.3) by —~ 
a 
du U.a , . Oa. a 
—— +—(u .V)u =-——Vp +vV'u (3.4) 
at Vv Vv 
Hereafter the superscript * will be omitted for nondimensional terms. Define the 


Reynolds number based on the diameter, Reg = 





Bes) , and take the curl of Eq. (3.4) to 
V 


be able to eliminate the pressure term and introduce vorticity, @=Vxwu_ where the 


vorticity vector, @, 1S 


W = 1,0),+leWtip We. 
Because of the absence of any azimuthal fluid motion about the axis of symmetry only the 
®-component of the vorticity survives and the transport equation for the ®-component of 


the vorticity is 


<0, -[Vx(uxa,i, }i — (v:-— 1 Jo, (3-5) 


= _ 
> Re, r’ sin’ 8 





Convert the velocity in Eq. (3.5) by defining the axisymmetric Stream function (y) 





Vex v 1 A od ; Azimuthal direction (3.6) 
rsin®@ 


In spherical coordinates, Eq. (3.6) is obtained as 


1 oy. 1 oy]. 
= A a ra a 
r -, sin 8 mI, sy ~_ . — 


U=U_1, tUgl, tUgly (3.8) 





where 


As noted earlier u, = 0 in this study. 


l I 
In the free stream y = A Ur’ sin’ @ or y= 5 r’ sin’ 6 in nondimensional form. 


By using the definition of the stream function it can be shown that 


aw 1 1 = I ¥)(_— | 
o=| nore sin@ 06 \sin@ 06 /\rsin® (3.9) 


T2 
Define an operator D” as 














er Yok tie mes S) 
(3.10) 


a —— 
or? —s r? =~ \ sin 8 00 


to obtain the azimuthal component of the vorticity; 


a 
=— Dp 
i rsin@ a 





Call [V X(u  X @pie)].ig in Eq. (3.5) the nonlinear term G(r,6,t). 


By carrying out the calculations in spherical coordinates it can be shown that 





| 1 |dcy,D’y) ~— 
V = ———___——. 2 
X(UXWolg) =ts| Ar, LL) +2D°wL(y) 
where 
Li =cos@ 
were en 
ae at Ou 


a(Z, 7257 —«C OZ OZ, eZ, OZ, 


a(x,y) Ox dy dy ox 


Combine all terms to get 


Re, 1 eee 


0 — = = 
5, (D’w) + 2 | ow +25 yLty) |=B*y 


(3.11) 


(3.12) 


(3.13) 


(3.14) 


Equation (3.14) is called the Vorticity-Stream function form of the Momentum 


equation in axisymmetric spherical coordinate system. For computational purposes, which 


will become clear later, define a modified stream function C(r,6,t), which is related to the 


usual Stokes stream function (y) by the relation; 
y =rCsin@ 
Also for any function f(r,@), define a new operator D* such that 
D?(fr sin®) =r sin@D7f 


where 


(3.15) 


l 





= - —— 
r’ sin? 9 
This gives 
o, =-D’C (3.16) 
u=Vx(Ci,) (367) 
u.= J <(Csin®) 
where rsin® 08 (3.18) 
0, =a 
: Tor 


The modified stream function C is written as the sum C=C+c, where 
— |. a: | - 
C= AS sin8, (which is the prescribed uniform-flow free stream condition), and c 
represents the disturbance produced by the presence of the rigid sphere. 


Re, 


Returning to Eq. (3.5), rearrange by multiplying both sides with r7 





> 





Re, do Re 
meee = G(r,9,t)+r’D’*o@ (3.19) 
r°D*C=-r’@ (3.20) 
where the nonlinear term 
G(r,6,t)=[Vx (ux @gig) ie (3.21) 
and the operator D7” is defined as 
pole 2 oO a : nics oes 3) 
~r? ort’ Or) fsin2000\ 00) r2sin?6 ae 


B. BOUNDARY CONDITIONS 


The flow is determined by numerically solving (3.5) for the vorticity and (3.20) for 


the modified stream function C subject to the boundary conditions; 
u=0 at hs (3523) 
u=U(t) “free stream velocity” as roo (3.24) 


Boundary conditions (3.23) and (3.24) may be reformulated in terms of C as 


dC 
ee, ee on r=] (3.25) 


] 
OQ=0;, C= > Ur sin(@) aS [— 0 (3.26) 


The Neumann Boundary Conditions in Equation (3.25) are handled with the 
Green’s function method that will be explained later. Homogenous boundary conditions 
on the axis of symmetry (8 = 0, 7) are automatically satisfied with the choice of the sine 
series expansions in the numerical method as shown later. 

C. BOUNDARY LAYER BEHAVIOR 

The boundary layer thickness on the surface of the sphere can be calculated with 

the method of power series expansion as discussed in, Schlichting (1987,pp.235-238). The 


body contour of the sphere with radius a is defined by r(x) =asin(x/a) where x is 


the distance along the surface from the stagnation point. In the boundary-layer limit the 
— 3 
velocity distribution near the surface of the sphere is given by U(x) = 5 U,, sin(x /a). 


By introducing the stream-function into the boundary-layer equations the final 


form of the transformation is obtained for axially symmetrical case; 











dy d° oy | 0° 0° 
oe (We SS es 3.27) 
dy dxdy \ox rdx /dy dx dy 
with the boundary conditions 
0 
ie y=0, S=0 (3.28a) 
dy 
0 
26 Sie (3.28b) 
oy 
The potential flow is defined by the series expansion as 
RO eee eX PUK +o csesscens Cotes (3.29) 
h ss aus “ fee). i y {3U a 
where = i = Sane ==> 
2 a aV v 
and the body contour is expanded to series with the same way; 
96.0 8 0 0nas Ob Gin e oo Oi ee (3.30) 


The stream-function is represented by the Blasius series in analogy with (3.29) & 


(3.30) as; 


W(x, y)= 5.7 {unaf, (n) + 2u5x"f, (n+ LSP er: } (3.31) 
| U, | 


By introducing equations (3.29),(3.30) and (3.31) into equation (3.27) we obtain a 


set of differential equations for the functions f, ,f,,....... 


If the first two terms of the expanded series are concerned the first equation is 
a | ae 
f, =-f,f, +o hy —1) (3.32) 
where differentiation with respect to n is denoted by primes. The boundary conditions are; 
1 Oe te (3.33a) 


10 


and the second equation is 


1 


i = Diet. = Pt ite ae aff = f,f; =|. (3.34) 
with boundary conditions 
n=0 f,=f, =0 (3.35a) 
ae 
=< f, = A (3.35b) 


After solving the Eq. (3.32) and Eq. (3.34) for f, and f,the stream function and 


the velocity can be calculated by using Eq. (3.31). In the boundary layer limit the 
boundary-layer thickness (at a given angular location) is chosen to be the radial distance 
from the surface of the sphere where the velocity reaches to its peak value. By using the 


definition of n it can be shown that 


[3 
Me. = Cen =) Deke (3.36) 


where r,.,, is the nondimensional distance where the velocity reaches its peak value. 


By using Eq. (3.36) the calculated 1,.,, values through the use of foregoing 
explanations are compared with the corresponding f,,, values which are obtained from 
the numerical results. The 1),,,, values are valid as Re, — co and they are compared 


with the boundary-layer thickness obtained from the numerical results. 








IV. NUMERICAL REPRESENTATION 


A. COMPUTATIONAL METHOD 

A numerical scheme called the pseudospectral or collocation method based on a 
recent paper by Chang & Maxey(1994) which is in term based on a technique developed 
by Marcus & Tuckerman (1987) is used to solve the Navier-Stokes equations. Spectral 
methods involve representing the solution to a problem as truncated series of known 
functions. From this point of view, spectral methods may be viewed as extreme 
development of the method of weighted residuals (MWR). The trial and the test functions 
are the key elements of the MWR. This 1s also true for the spectral methods. The key 
difference lies in the choice of the trial and test functions. The trial functions (also called 
the expansion or approximating functions) and the test functions (also known as weight 
functions) are chosen as infinitely differentiable global functions in spectral methods. 
Functions are transformed between physical space and spectral space through the use of 
Fast Fourier Transform (FFT) and Inverse Fast Fourier Transform (IFFT) which are 
numerically cheap and quick to compute. Although they are not exact due to aliasing and 
truncating errors they rapidly converges and few frequency modes are required to obtain 
accurate results. 

A physical process can be described either in time domain, by the values of same h 
as a function of t, e.g., h(t), or else in frequency domain, where the process is specified by 
giving its amplitude H, as a function of frequency f, that is H(f), with -co < f < o. If h(t) is 
real and even then H(f) is real and even or if h(t) is real and odd then H(f) is imaginary and 
odd. The processes can be transformed to either spectral domain or physical domain by 
using the Fourier transform equations. As it will be explained later in this chapter, we are 
dealing with odd functions by using sine series expansions for @, c. The properties of odd 
and even functions are also valid for their derivatives and these properties will help to 


increase the computational efficiency in the calculations. 
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The functions which are introduced into the governing equations are presented 
both in spectral space as a finite series of trial functions and by values at collocation grid 
points in physical space. The vorticity @ and the modified stream function C are expanded 
as a Chebyshev polynomial in the radial direction and as sine series in the azimuthal 
direction. The calculations in the radial direction are done in physical space while the 
derivatives in azimuthal direction are obtained in spectral space. Since each sine term in 
the expansion satisfies the homogenous 8 boundary conditions (@ = 0, ‘Y = 0 at @ = 0,7) 
exactly, no further @ boundary conditions need to be applied. 


In general let a variable, say @, be expressed as 


N 
o(r,6,t)= > @, (z,t)SinO, (4.1a) 
with 
@,(z,t)= 4, (t)T,,(Z) (4.1b) 
m=0 


where T,,(Z) is Chebyshev Polynomial with -1 <z <1 





9 ee 
Z ) 


An algebraic map is used to map the radial interval 1 <r <r.. to the interval 


-1 <z< 1 where 








] 
1+z \e 
r= [i ay pe izI<1 (4.4) 
b-z 
2 
with b=l+— (note that r, = 1, the surface of the sphere) 
r-f, 


L is a scaling parameter which is used to map the radial interval to the z-interval. A 
finite, but large outer r.. was chosen to avoid regularity problems in the radial 


differentiation. Another mapping parameter,a is defined for higher Reynolds numbers 
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(Re, 220). For small Reynolds numbers @ is taken as 1. Two typical grids for two 


different parameters are shown in Figure 1. 
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Figure 1. Typical grids for @=1 a)L=8 rinf=250 M=64 N=64 b) L=2 rinf=50 M=64 N=64 


The sine expansions in Eq.(4.1a) and Eq.(4.1b) for @ and C satisfy the 
homogenous 8 boundary conditions and match exactly the periodicity and symmetry 
conditions in 8. The solution of the elliptic Eq.(3.20) is required to obtain the modified 
stream function C from vorticity @ at any time level. Following Marcus & Tuckerman 


(1987) separable derivative operators D,’ and De” are defined as 


re — +— Di (4.5) 
) ) 
lbs = 2 il 
"Or (1 =| 
where (4.6) 


D.= sino 5 (sino 5) - 1 


In spectral 6-space the effect of operator De’ at any fixed radial location (physical r-space) 


can be written as 





D,.f£(8) -| sin 0 2 “+ sin OO — VE, sin(j8) 
00? 00 
(4.7) 
= y J - (cba > yj’) sin( j®) + ; j(1 + j) sin(28 + 38) + id — 7) sin(20 — | 
j=! 


If we write this expression in the form of Déf(@) = /S i ii ; we get 


wes PPE Gece: 1=j+2 
i” . 
S,= oe ee 1=] (4.8) 
G+2)0+1) 
A tetereseeeees 1= j-2 
Oe Berane sc... ee otherwise 


This matrix is NxN if the representation for D¢f is truncated to N Fourier terms. 


Similarly the operation of multiplying a function f(®) by sin“ @ produces another NxN 
matrix which is called A matrix where the only non-zero terms are ; 


ag ia HAIG + Doses 53 


— (4.9) 
ACC 7 1< j,1+ j=even 


B. TIME INTEGRATION 

Time integration of the Eq.(3.19) and Eq.(3.20) are accomplished through the use 
of an explicit second-order Adams-Bashforth scheme for the nonlinear terms and an 
implicit second-order Crank-Nicholson scheme for the viscous and linear terms. The 
calculations are made in physical r-space and spectral 0-space. If we denote vorticity a@r,t) 


as the vector of N sine series then discretized forms of Eq.(3.19) and Eq.(3.20) are; 


Bes] Soa. , Re, ; (S So 
ro | —— [=r —“]0G, + BG, |+|D, + Aj} ———_— 4.10 
’) At a t B sol | 2 ( ) 


Rearrange Eq.(4.10) and get ; 





R R 
b: Aer vi ow =: +A¢r? Ret bo ~r? Re,[aG, -BG,_,,] (4.11) 


3 l 
where a=73B=-> and G(r,0,t)=[Vx(ux@)]i, 


Equation (3.20) is also discretized by using the same method as ; 
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[D? + Ale(r,t+ At) =—r?@(r, t+ At) (4.12) 


The radial derivatives are evaluated by collocation methods and expressed as 
matrix operations on the vector of function values at the grid points in physical r-space. 


For (M+1) collocation points, we introduce the (M+1)x (M+1) matrix operation, that is 
R 
Ey Gye D- +a -Reee} (4.13) 
t 


where, I is the identity matrix. The operator [D,’-r°(Re ,/At)+A] in Eq.(4.11) may now be 


written in block matrix form as 


Da o@8 aor 4 
0 D®2) oO a22)I 
0 0 DB) 0 
0 0 0 p®(4) 


(4.14) 


This matrix acts on the vector of a length of (M+1) x N for vorticity coefficients. 
Fq.(4.11) and Eq.(4.12) are upper triangular, block matrix problems in physical r-space 
and spectral @-space that can be solved by any solver with the Dirichlet boundary 
conditions discussed below. The Eq.(3.32) is treated in the same manner and leads to an 
upper triangular, block matrix problem with each diagonal term on row n replaced by 
D,’+al(n)I. 

C. NONLINEAR TERMS 


If one opens the nonlinear term G(r,0,t) in spherical coordinates gets six terms; 

















_ 190, ac 10®, ac _cot(®) ac 
r or oO r 00 or r “aie 
(4.15) 
_ 20t(8) -, 9% ®, AC C 00, 
r Or r’ 00 r’ 00 


The nonlinear terms in Eq.(4.15) are handled by transforming the spectral 


coefficients to physical space first and transforming the result of the products to spectral 
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space after the multiplication. Basically the derivatives are obtained in spectral space and 
the products are done in physical space and finally the result is transformed to spectral 
space. 
D. GREEN’S FUNCTION METHOD 

In order to enforce the radial boundary conditions in Eq.(4.11) and Eq.(4.12), 
Green’s function method has been developed which gives the correct boundary conditions 
in Eq.(3.25) and Eq.(3.26) and allows the surface vorticity to develop naturally. If @ and 
c consist of two parts, homogenous and particular solutions, we write 
o=0,+>26, and c=, + YL j= ee N (4.16) 

j=l j=l 

If we introduce Eq.(4.16) into Eq.(4.11) and Eq.(4.12), we find the homogenous parts 


with corresponding boundary conditions as follows; 


E*® (1,6) =0 (4.17) H’t,=-r’@, (4.18) 
@ (r =1,0;) =9, C,(r =1,0;) =0 
0 (r=,0;) =0 E(r=,6;) =0 


The particular parts with corresponding boundary conditions are ; 


Eros (r,8,t + At) = R(r,9,t) (4.19) Hoon Ga! t+At)= —r°@, (r,9,t) (4.20) 


® ,(r =1,0) =0 c,(r=1,8) =-Cl,_, 
@ , (r= 0,8) =0 Cc (rt — cong — 0 
P 
E’ and H’ are obtained from block matrix defined as in Eq.(4.14) 
aC 
To find A, we enforce a =——l_, (4.21) 
or Or 
If we substitute Eq.(4.16) into Eq.(4.21), we get 
dc NOC. dc 
=— ape | 2) see | (4.22) 


= 
or r=1,0; j=! Or r=1,0; or r=1,6; 





c oc dC. 
IES = e a =H, and ——~ | =A, (4.23) 
or Ola or nn 


N 
finally from Eq.(4.21),Eq.(4.22) and Eq.(4.23), we find H,= >) A,A, and solve for 


jel 
rN ; as X = A"H. These A values can now be used in Fq.(4.16) to arrive at the values of w 


and c. 
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V. DISCUSSION OF THE RESULTS 


The results are going to be discussed in three parts. The drag coefficients 
calculated by following the explanation below for various Reynolds numbers are compared 
with the results of Chang & Maxey (1994) and the streamlines contours are compared 
with the experimental flow visualization pictures and finally the boundary-layer behavior in 
the boundary-layer limit is discussed by using the numerical results. 

A. EVALUATION OF THE DRAG 


By integrating the contributions from the stress tensor 0, over the sphere surface 

the resultant fluid force on the sphere exerted by the moving fluid is found as; 
F, = $o,n,dS (5.1) 
where n is the unit outward normal. There is no lift force because the flow is axisymmetric 


and the force F acts parallel to the free-stream velocity U. The components of the force 


can be expressed in terms of the spherical polar coordinates as ; 


F, = $(6,, cos(0)— 6 sin(®))dS (5.2) 
The normal component of the stress 6, is 


2 pee (5.3) 
0, =—p+2pv—> | 


The no slip boundary conditions in Eq.(3.23) and the condition of incompressible flow in 


Eq.(3.2) causes the du,/dr to vanish on the surface of the sphere and the only 


contribution from o_ is through pressure. The shear stress 6, is 
tao) (ue 1 ou, 
= Se | > | 5.4 
Ore 25 2(*4] 2r a oo 


/| 


Equation (5.3) and Eq.(5.4) are in dimensional form. By using the non-dimensional 


variables the drag force C, is defined as 


C,=byapa U- (5.5) 
The drag force C,may be split into two separate parts, C,for the frictional 


component due to shear stress and C,for the pressure component. The frictional 


component is calculated from Eq.(5.2) and Eq.(5.4) as 
€ walla 0) sin* (8)d0 (5.6) 
: Rey 


The pressure component C, in non-dimensional form 


C, =—| p(r = 1,6)sin(26)d6 (5.7) 
0 


The pressure variation over the surface of the sphere can be determined once the 


vorticity field is known : 


270 
p(r = 1,8) = p(1,m)— = J = (ro) | a8 (5.8) 
0 r=1 


A non-dimensional pressure coefficient k(@) can be defined to make the 
calculations easier; 

k(9) = 2p(r = 1,8) (5.9) 

The pressure at the stagnation point, @=7,can be evaluated from the radial 


momentum equation as; 


k nmacein aden | dr (5.10) 
Re +r 00 e=n 
B. STEADY FLOWS 
In order to verify the numerical method and to provide a comparison, simulations 
for steady, unidirectional flow were carried out for Re<104. The parameters which 


were used in the calculations are listed in Table 1. 


Ze 


For all calculations the flow is initially at rest and then a uniform, constant external 
flow is introduced in the direction of @=7 and the simulation continued till the solution 
has converged to steady state. Two different grids were used in the calculations. A spatial 
discretization of 64x64 grid points was used for Re <10 and a relatively finer mesh of 
128x128 grid points was used for Re2 20. For Re220 a mapping parameter “a” is 
defined for finer resolution near the surface. The vorticity diffuses over a large region at 
low Reynolds numbers compared to the sphere radius while it is convected downstream 


with a more complex structure, which needs increased resolution. 


a a 
oe or 
_ 0.05 | 250 
a 
os fg of 005250 
pos 0.05) all 250 
00s > a0 
“oe te te oe 
09 8 pet 0.05 250 

. a 
| 002 | 50 
= es a 
meee | 158 is) i ee 


Table 1. Parameters used in the calculations 
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The present results have been summarized in Table 2 with the previous results 
which were calculated by Chang & Maxey (1994). As it can be seen in Table 2 there is a 
general agreement between the present results and the previous results. The separation 


angle,@., and the ratio of bubble length (s) to the diameter of sphere (d) are plotted later 


in this chapter to compare them with the results which were obtained experimentally by 


Taneda (1956). 
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Figure 2 shows the variation of vorticity and streamlines in the vicinity of the 


sphere at Re = 0.1. As it can be noticed the vorticity is diffusing over a large region and 


very slight asymmetry about 8 = 1/2 is observable. 





Figure 2. (a) Lines of constant vorticity for Re = 0.1 A@ = 0.1 (b) Streamlines Ay = 0.3 
Figure 3 shows the constant vorticity and streamlines contours for Re = 5 .The 


asymmetric region 1s getting more observable about 86=7/2 as Reynolds number 


increases. 


4 





Figure 3. (a) Lines of constant vorticity for Re = 5 A@ = 0.15 (b) Streamlines Ay = 0.25 


In Figure 4 the constant vorticity and streamlines contours are plotted for Re = 40. 
A separation region is clearly observable on both plots. The recirculation in the separation 


region is also noticeable with the asymmetry over the sphere. 


ye, 





0 
—2 0 2 4 
Figure 4. (a) Lines of constant vorticity for Re = 40 Aw = 0.1 (b) Streamlines Ay = 0.05 


for ‘----‘, Aw = 0.005 for ‘-.-.-.-‘, Aw = 0.00004 for the separation region. 


Figure 5 magnifies the separation region over the constant streamline contours for 
Re = 40. There is a distortion of the vorticity field caused by advection in the flow and a 


separation is clearly visible with a small region of negative streamlines in recirculation. 





—2 -1.5 —1 -0.5 0 0.5 1 1.5 2 25 3 


Figure 5. Constant streamline contours for Re = 40. Ay = 0.00004 in the separation 
region , A = 0.005 for ‘-.-.-.-‘, Ay = 0.2 for ‘- - - - ‘ 
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Figure 6 shows the constant vorticity and streamlines for Re = 56.5. The 
separation and the recirculation inside the separation region is clearly visible. There is a 
noticeable change in the wake length as the Reynolds number is increased from Re = 40 to 


Re = 56.5. 
4 F 





Figure 6.(a) Constant vorticity contours for Re = 56.5 A@ = 0.3 (b) Streamline contours 
Aw =0.005 for ‘----‘, Aw = 0.15 for ‘-.-.-.-‘, Aw = 0.0003 for the separation region. 
The constant vorticity and streamline contours for Re = 104 are plotted in Figure 
7. The separation region and recirculation inside the separated flow is clearly visible. The 
increase in the wake length is also noticeable for this particular Reynolds number. It is also 
noticeable from the foregoing plots that the vorticity diffuses to smaller region as 


Reynolds number increases. 


4r 





Figure 7. (a) Lines of constant vorticity for Re = 104 A@ = 0.4 (b) Streamlines Aw =0.005 
for ‘----‘, Ay = 0.15 for ‘-.-.-.-*, Aw = 0.0009 for the separation region for Re = 104. 
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Figure 8 magnifies the streamlines for Re = 56.5. The flow visualization picture 
for same Reynolds number is also presented in Figure 9. It is clearly observable that there 
is a simuarity between Figure 8 and Figure 9. The wake length, separation location and the 
nature of the streamlines in Figure 8 and Figure 9 look alike in both figures for the 
numerical results and flow visualization picture. 

Figure 10 magnifies the streamlines plot for Re = 104 to present the similarity of it 
with the flow visualization picture in Figure 11 for the same Reynolds number. The 
extension of the recirculating wake to a full diameter in the downstream 1s clearly visible in 
Figure 10 and Figure 11. The flow is still laminar at this Reynolds number. 

The calculated separation angles and the ratio of the wake lengths to the diameter 
of the sphere are presented in Figure 12 respectively along with the experimental data 
which were obtained by Taneda (1956). There is a good agreement between the present 
results and Taneda’s data. Taneda has also extrapolated his data and conjectured the 
critical Reynolds number of the “onset of separation” as Re, = 24 .On the other hand, the 
earlier work of Nisi&Porter (1926),which was specially designed to detect the onset of 
separation with ultramicroscopic methods, reported separation at Re = 10 and indicated 
that the onset of separation in an infinite fluid is about Re, =8. 

An estimate of Re = 20 is given by the present results as the Reynolds number for 
the onset of separation. This is good agreement with the value of Re, = 20.7 given by 
Chang&Maxey (1994), the value of Re, = 20.5 given by Dennis& Walker (1971) and the 
value of Re, = 20 given by both Le Clair et al.(1970) and Lin&Lee (1973). Others have 
found separation to first occur at different Reynolds numbers, for example Rimoné&Cheng 
(1969) found separation as low as Re, = 10. 


Taneda’s data also suggests that the length of the separation bubble varies linearly 
with the logarithm of the Reynolds number. The present results are consistent with such a 
result and they remark that over this limited range a simple linear relationship with 


Reynolds number is a good approximation. 
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separation region for Re = 56.5 
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7 Figure 9. The flow visualization picture of streamlines for Re = 56.5 ( Archives de 
1” Academie des Sciences de Paris. Payard & Coutanceau 1974 ) 
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-2 —1 0 1 | 2 3 | 4 
Figure 10. Streamlines Ay =0.005 for ‘----‘, Aw = 0.15 for ‘-.-.-.-‘, Aw = 0.0009 for the 
separation region for Re = 104. 








Figure 11. Visualization picture for Re = 104 by a thin coating of condensed milk on the 
sphere, which gradually melts and is carried into the stream of water. Taneda 1956. 
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Figure 12.(a) Separation angle 0. vs. Log Re (b) The ratio of bubble length to diameter 
(s/d) vs. Log Re ‘---* Taneda (1956) 


C. BOUNDARY LAYER BEHAVIOR 
Once the stream function was known for the particular flow the velocity over the 


sphere can be calculated by using the Eq.(3.7). The velocity in @ direction u, which is 


parallel to the surface of the sphere can be calculated by using the relation 





1 oOo 
ig (5.11) 
rsin@ or 
We also know from potential flow theory that the stream function is given by 
1 leslie 
vas(r 2) sin’ Q (5.12) 
Z i 


By using the relation in Eq.(5.11) u,can be calculated from Eq.(5.12) for potential flow 


as 





&) 
u, = (2 Jain (5.13) 


Figure 13 shows the u,/sin® values which are obtained from Eq.(5.13) for 


potential flow along with the numerical results for various Reynolds numbers from this 


study. It is clearly observable that the velocity profile matches with the potential flow at 
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higher Reynolds numbers. The velocity profile for a small Reynolds number is also 
presented in Figure 13 for comparison, and shows that the velocity reaches its free stream 
value over a much larger distance in this Stokes flow regime. 


1.5 


Free stream value 


Ug /sin8 


0.5 





; 
1 1.05 1.1 1.15 1.2 1.25 1.3 1.35 1.4 1.45 7.5 
I 


Figure laitie: ue / sin9 profiles for varying Reynolds numbers (8 = 1 / 4from the nose) 
‘-----* represents potential flow, solid lines represent numerical results. 


The approach to boundary-layer behavior is now examined by comparing the 


boundary-layer thickness from numerical results with values of 1),.,, from boundary-layer 
analysis as discussed earlier in Chapter II.C. Figure 14 shows the 1. values at 
corresponding Reynolds numbers along with the 1),.,, value as Re— ce for a particular 
angle. 


The 1... Values for corresponding Reynolds numbers can be modeled by assuming 
a model such as 





C 
=C,+— 5.14 
Saar (5.14) 


where C, is the 1],.,, value as Reo ,C,, and m are the constants. 
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10 


Log(Re,) 
Figure 14. 1),.., values vs. Log Re ‘---" represent the 1],.., value as Reco 
By using the ics values at Re=30, Re = 56.5 and Re = 104 , C, and mcan be 
calculated from Eq.(5.14) as C, = 2.464 and m = 0.182 to obtain the model as 


2.464 


Reo (5.15) 


Q=C,+ 
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VI. CONCLUSIONS AND RECOMMENDATIONS 


The numerical results have showed a good agreement with the other results which 


were obtained by using different numerical methods. The pressure coefficient, C,, which 
is one of two components of the total drag coefficient, C,, couldn’t be resolved for higher 


Reynolds numbers through the use of present method. But the calculated skin friction 


coefficients, C, , and the stagnation-point pressure k(7) have been resolved successfully. 


A fine mapping for collocation points is very important to obtain more accurate 


results. The mapping parameters r_,Q@ and L should be chosen in a way to provide 


enough collocation points in the boundary layer and recirculation regions. A finer mapping 
can increase the accuracy by allowing finer resolution in these large gradient regions. 

The Green’s function method that provides the natural development of vorticity 
along with the Neumann boundary conditions should not be considered as a trivial detail 
and must be used to have a stable convergence process. 

A different FFT technique can be used which does not restrict the number of grid 
points to a power of two. This could reduce the memory requirement by allowing a 
relaxation in the choice of grid points. 

Very good agreement has been obtained with the previous calculations and 
available experimental flow visualization pictures. This study is at a stage where it can 


now be extended to explore higher Reynolds numbers with suitable turbulent models. 
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APPENDIX A. PROGRAM STRUCTURE 


In the first part of the program, problem parameters were defined such that Re, dt, 
L, r_inf, M, N. Calculations which don’t depend on time were also performed in this part. 
In order to gain some additional memory D1, D2 and D3 matrices in (4.11)&(4.12) were 
calculated with the function D_Bmat and saved in terms of diagonal elements only. Non- 
diagonal elements of these matrices were included separately in each iteration. 

Calculations showed that the matrix M1 (derivative operator in r direction) and 
consequently the matrix M3 and the block matrices which are formed by the diagonals of 
D1, D2 and D3 are ill-conditioned. To avoid the numerical errors in solving the systems of 
equations, Singular Value Decomposition (SVD) method was used. SVD of the matrices 
and calculations of w and c that are used in Green Function method were performed in this 
part. The built-in function SVD in MATLAB, was used for decompositions of the 
matrices in the systems of equations. 


The calculation part basically consists of four steps. In first step, A,; matrix was 


calculated from eqn. (4.22) out of the time loop. In second step, inside the time loop the 
particular solutions of w and c were found from eq. (4.17) and (4.18). In third step, H that 
is used in Green’s Function method, was calculated and finally homogenous and particular 
solutions are combined to get w and c. 

Most of the calculations were performed in spectral domain. Only nonlinear 
Operations were done in physical domain. Nonlinear terms were calculated with the 
functions nonlin for w-c equations . 

A convergence criterion isn’t used in the program to avoid the wrong results which 
are caused by the over-shooting for high Reynolds numbers. Instead, the results are 
checked during the calculations and program was terminated after seeing the steady state. 

Solsvd : This function is used to solve equations after SVD in Matlab. A limit, 
which is defined as a number smaller than the values of S matrix in SVD, is used in this 


function to avoid the numerical errors due to those ill-conditioned matrices. This is 


o/ 


avoided by forcing those singular values which are smaller than the defined limit, to 
zero. 10°'* is used in all calculations. 

Solsvdl : Since this function was used for solving homogenous parts in Green’s 
Function method and performed once before the time-loop, no limit was used in solsvd/. 

Phy : This function was used for transformations of the variables from spectral 
domain to physical domain 

Operat : Derivatives in 8 direction in spectral domain was performed with the 
function, operat. Since these terms were used only in nonlinear terms the output from this 
function is in physical domain. 

L __s: Mapping parameter 

Re —_: Reynold’s Number 

dt : Time interval between iterations 

r_inf : Outer boundary 

step : Number of iterations 

N — : Number of grid points in @ direction. 

M — : Number of grid points in r direction. 

xcor, ycor : Cartesian coordinates of grid points with respect to center of sphere. 
(used for presentation of the results) 

M1 _ : Derivative operator in r direction. 

D1, D2, D3 : Matrixes defined as in (4.11)&(4.12) 

0) : Vorticity. 

€ : Potential function as in (3.15a)&(3.26) 

wtil, ctil : Homogenous part of the w and c defined as in Green Function. 

Aij, H, landa : Variables defined as in Green’s function method in chapter IV. 


wp, cp: Particular solution of @andc defined as in Green Function. 


Cf — : The frictional component of drag 
kpi § : Non-dimensional pressure coefficient at 6 = 7 
ep : Pressure component of drag 
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APPENDIX B. PROGRAM CODES 


%% define program parameters 

% Table 1 list these parameters for every Reynolds numbers 
clear all 

global E M1 M2 M3 

L=8;Re=1;dt=0.02; 

r_inf=250;step=10000; 


%7% define grid points 


N=32;n=[0:N-1]’; 
teta=n.*pi./N; 
M=32;m=[0:M]’; 
Z=cos(pi.*m./M); 
b=1+2*L/(r_inf-1); % map=1 
r=1+L.*(1+z)./(b-z);% map=1 
global NM 

xcor=cos(teta)*r’; 
ycor=sin(teta)*r'; 


E=emat(M); 
disp(‘E matrix is done’) 


%% define M1,M2,M3 


dz_dr=(b.*(L+r-1)-b.*(r-1)+L)./(L+r-1).%2; % map=1 
B=diag(dz_dr); 

M1=B*E; 

R=diag(r);global R Re dt 

M2=R*M1; 

M3=R*(2.*eye(M+1)+M2)*M1; 


%% define D1 D2 D3 
[D1,D2,D3]=D_Bmat(M,N); 
% define initial wc wp cpu 
dum1=(-1.5./r.42)*sin(teta’))’; 


dum2=[dum1;zeros(1,M+1);-1.* flipud(dum1(2:N,:))]; 
dum3=fft(dum2); 
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w=imag(dum3(1:N, :))];w=w’; 

dum 1=(-0.75+0.25./r.*2)*sin(teta’))’; 
dum2=[dum1;zeros(1,M+1);-1.*flipud(dum1(2:N,:))]; 
dum3=fft(dum2);c=imag(dum3(1:N,: ))];c=c’; 
wp=zeros(M+1,N);cp=wp; 


% define kmat (for teta der.) and cotteta (for nonlinear terms) 


dum=ones(1,M+1); 
kmat=n*dum; 
cotteta=[0;cot(teta(2:N))]*dum; 
global kmat cotteta teta 


% svd decomposition for D1 D3 and define rsquare (nonlinear) 


D lu=zeros(N*(M+1),M+1); 
Dil=D1u;D1p=D1u;D31=D1u;D3u=D1u;D3p=D 1u; 
rsquare=zeros(M+1,N); 
for 1=1:N 
bas=(1-1)*(M+1)+1; 
son=bas+M;; 
rsquare(:,1)=r.‘2; 
D1(bas,:)=[1 zeros(1,M)]; % modify D1 for w @ r=inf 
D1(son,:)=[zeros(1,M) 1]; % modify D1 for w @ r=] 
[D11(bas:son,:),D lu(bas:son,:),D 1p(bas:son,:)]=svd(D 1(bas:son,:)); 
D3(bas,:)=[1 zeros(1,M)]; % modify D3 for c @ r=inf 
D3(son,:)=[zeros(1,M) 1]; % modify D3 for c @ r=1 
[D31(bas:son,:),D3u(bas:son,:),D3p(bas:son,:)]=svd(D3(bas:son,:)); 
end 
global rsquare 


% define c_c_tetac_ru_u_tetau_r 


c_=(0.5.*sin(teta))*r’ ; 
c_teta=0.5.*cos(teta))*r’; 
c_r=-(0.5*sin(teta))*ones(1,M+1); 
global c_c_teta c_r 


% define c boundary at r=1 


dum4=-0.5.*sin(teta); 
dum5=[dum4;0;-1.* flipud(dum4(2:N,:))]; 
dum6=fft(dum5S); 
cbr1=imag(dum6(1,N))’; 
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% begin green function wtil ctil 


wtil=zeros(M+1,N‘2); 
ctil=wtil; 


for jj=1:N 

RR=zeros(M+1,N); 

dum 1=[zeros(1,jj-1) 1 zeros(1,N-jj)]; 

dum2=[dum1 0 -1.*fliplr(dum1(2:N))]; 

dum3=fft(dum2); 

dum4=imag(dum3(1:N));if jj==1;dum4=real(dum3(1:N));end; 

RR(M-+1,:)=dum4; 

wtil(:,jj*N)=solsvd1(D1]((N-1)*(M+1)4+1:(N-1)*(M+1)+14M,:).... 
D1lu((N-1)*(M+1)+1:(N-1)*(M+1)+14M,:).... 
D1p((N-1)*(M+1)+1:(N-1)*(M+1)4+14M,:),... 
RR(:,N)); 


for 1=N-1:-1:1 
sum=zeros(M+1,1); 
11=(jj-1)*N-+1; 
for j=i+1:N 
if rem(i+j,2)==2 
sum=sum+(-2*(i-1)).*wtil(:,(jj-1)*N-+j); 
end 
sum(1)=0; 
sum(M+1)=0; 
end 
RR_=RR(:,1)-sum; 


wtil(:,11)=solsvd1(D11((i-1)*(M+1)+1:-1)*(M+1)+14M,;).... 
D1u((i-1)*(M+1)+1:G-1)*(M+1)+14M,:).... 
D1p(G-1)*(M+1)+1:G-1)*(M+1)+14M,:),... 
RR_); 


end 


tempwtil=(-1.*rsquare).* wtil(:,Qj-1)*N+1:jj*N); 

tempwtil(M+1,:)=zeros(1,N); 

tempwtil(1,:)=zeros(1,N); 

ctil(:,jj*N)=solsvd1(D3]((N-1)*(M+1)+1:(N-1)*(M+1)+14+M,°),... 
D3u((N-1)*(M+1)+1:(N-1)*(M+1)+14M,:),... 
D3p((N-1)*(M+1)+1:(N-1)*(M+1)+14+M,:).... 
tempwtil(:,N)); 
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for i=N-1:-1:1 
N=(jJ-1)*N-+1; 
sum=zeros(M-+1, 1); 
for j=1+1:N 
if rem(i+j,2)==2 
sum=sum+(-2*(i-1)).*ctil(:,Qj-1)*N+)); 
end 
sum(1)=0; 
sum(M+1)=0; 
end 
wtil_=tempwtil(:,1)-sum; 


ctil(:,11)=solsvd 1(D31((i-1)* (M+1)+1:G-1)*(M+1)+14+M,;).... 
D3u((i-1)*(M+1)+1:(-1)*(M+1)+14+M,:).... 
D3p((i- 1)*(M+1)+1:G-1)*(M+1)+14+M,:),... 
wtil_); 
end 
end 


% find Aij matrix 

for jj=1:N 
blok=ctil(:,Gj-1)*N+1:jj*N); 
dum 1=phy(blok’)’; 
dum2=(M1*dum1)’; 
dum3=[dum2; zeros(1,M+1) ; -1.*flipud(dum2(2:N,:))]; 
dum4=fft(dum3);dum5=imag(dum4(1:N,:));dum5=dum5'; 
ctilr_1Qqj, 1:N)=dum5(M-+1,:); 

end 

Ay=ctilr_1’; 

% find svd of Ajj for landa solving 


Aij(1,:)=[1 zeros(1,N-1)]; 
[Aiju, Aijs,Aijv]=svd(Aij); 


% begin time loop 
for ops=1:step 


if ops==1 ;prc=c;prw=w;pru=u;end 
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% find Rl 


R1=zeros(M+1,N); 
R1(¢:,N)=D2((N-1)*(M+1)+1:(N-1)*(M+1)+1+M,:)*w(:,N); 
for 1=N-1:-1:1 
summ=zeros(M+1,1); 
for j=i+1:N 
if rem(i+j,2)==2 
summ=summ+(-2*(i-1)).*w(:,j); 
end 
end 
R1¢.,)=D2(i-1)*(M+1)+1:G-1)*(M+1)+14+M,:)*w(:,i)+summ; 
end 
Rer—=)..* Ri; 


% find R2 


R2a=nonlin(w,c); 
R2b=nonlin(prw, prc); 
R2=3/2.*R2a-1/2.*R2b; 
=(-1*Re).*rsquare.* R2; 
Rtot=R1+R2; 
prw=w;prc=c; 


%% green function solution 
% find new wp 
%% apply BC for wp 


Rtot(1,:)=zeros(1,N); % BC for wp=0 at r_inf 
Rtot(M+1,:)=zeros(1,N); % BC for wp=0 at r_1 


wp(:,N)=solsvd(D11((N-1)*(M+1)+1:(N-1)*(M+1)+14+M,:).... 
D1u((N-1)*(M+1)+1:(N-1)*(M+1)+14M,:),... 
D1 p((N-1)*(M+1)+1:(N-1)*(M+1)+14M,:).... 
Rtot(:,N)); 


for 1=N-1:-1:1 
sum=zeros(M+1,1); 
for j=i+1:N 
if rem(i+j,2)==2 
sum=sum+(-2*(i-1)).*wp(:,j); 
end 
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sum(1)=0; 
sum(M-+1)=0; 

end 
Rtot_=Rtot(:,1)-sum; 


wp(:,1)=solsvd(D11(G- 1)*(M+1)+1:G-1)*(M+1)+1+M,:).... 
D1u((i-1)*(M+1)+1:G-1)*(M+1)+14M,:).... 
D1p(G-1)*(M+1)+1:(-1)*(M+1)+14+M,:).... 
Rtot_); 


end 
% solve for new cp 


tempwp=-1.*rsquare.* wp; 
tempwp(M+1,:)=cbr1; % BC for cp=-c_r at r_1 zero for pure spin 
tempwp(1,:)=zeros(1,N); % BC for cp=0 at r_inf 


cp(:,N)=solsvd(D31((N-1)*(M+1)+1:(N-1)*(M+1)+14+M,:).... 
D3u((N-1)*(M+1)+1:(N-1)*(M+1)+1+M,;).... 
D3p((N-1)*(M+1)+1:(N-1)*(M+1)+14+M,:),... 
tempwp(:,N)); 


for 1=N-1:-1:1 

sum=zeros(M+1,1); 

for j=i+1:N 
if rem(i+j,2)==2 

sum=sum+(-2*(1-1)).*cep(:,J); 

end 

sum(1)=0; 

sum(M+1)=0; 

end 

w_=tempwp(:,1)-sum; 


cp(:,1)=solsvd(D31((i- 1)*(M+1)+1:(-1)*(M+1)+14+M,:).... 
D3u(d-1)*(M+1)+1:G-1)*(M+1)+14+M,:).... 
D3p((i-1)*(M+1)+1:(1-1)*(M+1)+1+M,:).... 
w_); 
end 


% find landa 
% find H 


dum1l=phy(cp’)’; 
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dum2=(M1*dum1)'; % NxM+1 

dum3=-1.*c_r-dum?; 
dum4=[dum3;zeros(1,M+1);-1.*flipud(dum3(2:N,:))]; 
dum5=fft(dum4); 

dum6=1mag(dum5(1:N,:)); 

dum6=dum6'; % M+1xN 


=dum6(M-+1,:);H=H'; % Nx1 
Fi(1)=0; 
landa=solsvd(Aiju,Atjs,Aijv,H); 


% find complete w & c 


dum1|c=zeros(M+1,N);dum1w=dum1c; 

for jj=1:N 
blokc=ctil(:,Qj-1)*N+1:3j*N); 
blokw=wtil(:,Qj-1)*N+1:j)*N); 
dum2c=landa(jj).*blokc; 
dum2w=landa(jj).*blokw; 
dum1c=dum1c+dum2c; 
dum1w=dum1w+dum2w; 

end 


W=wp+dum1w; 
c=cp+dumlc; 


res=phy(w’); 
resp=phy(prw’ ); 
resc=phy(c’); 
resc=resc+c_; 


% Calculate Cf for every iteration 
dum 1=res(:,M+1); 
dum2=(sin(teta)).%2; 
dum3=dum1.*dum2; 
cf(ops)=-4/Re*trapz(teta,dum3); 


% Calculate k(7) for every iteration 


dum 1l=w’; 
for 1=1:M+1 
sum(1)=0; 
for j=1:N 
f=((j-1)*duml (i,j).*(-1.4G-1)))./N; 
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sum(i)=sum(i)+f; 
end 
end 
dum2=(1 ./r).*sum; 
kpi(ops)=(1+(8/Re)*trapz(r,dum2)); 


% Check for convergence 


if max(max(isnan(w)))==1;error(' it blowed up!!!! ');end 
% Save the variables to a file 
save tez.mat c w Re ops r_inf cf kpi LN M dt; 
end % end of time loop 
Jo To [o To Vo 


function [D1,D2,D3]=D_Bmat(M,N) 
global M3 R dt Re 
k=(M+1)*N; 
D1=zeros(k,M+1); 
D2=zeros(k,M+1); 
D3=zeros(k,M+1); 
for 1=1:N 
for j=1:N 
0 
D1(G-1)*(M+1)+1:4-1)*(M+1)+14M,:)... 
=M3+(-1*(1-1)*(1).*eye(M+1)-Re/dt.*R.%2); 
D2(-1)*(M+1)+1 :(7-1)*(M+1)+14+M,:)... 
=M3+(-1*(1-1)*(1).*eye(M+1)+Re/dt.*R.%2); 
D3((i-1)*(M+1)+1:(1-1)*(M+1)+1+M,:)... 
=M3+(-1*(i-1)*(1).*eye(M+1)); 
end 
end 
end 
Jo To To To To Yo Yo 
function E=emat(M) 
E=zeros(M,M); 
m=0:M; 
Z=cos(pi.*m./M); 
for 1=0:M 
for j=0:M 
if 1~>] 
if (i==0 | 1==M) & (G==0! j==M)) | (G~=0 & i~=M) & (j~=0 & j~=M)) 
c=]:;else 
if (G==0 | i==M) & (j~=0 & j~=M) 
c=2;else 
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if (i~=0 & 1~=M) & G==0 | j==M) 
c=1/2;end;end;end; 
E(it1,j+1)=c*(-1)4Ci+j))(z+1)-zG+1)); 
else; 
it i=—0 
E+] j+1)=(2*M‘2+1)/6; 
else 
if i==M 
EG+1,j+1)=(2*M‘%2+1)/(-6); 
else 
EG+1,j+1)=-1*z9G4+1)/(2*(1-z(j+1)*2)); 
end 
end 
end 
end 
end 
JoTo To To To To Go 


function R2=nonlin(w,c) 

global kmat rsquare cotteta N M M1 c_c_teta c_r teta 

% variable list 

Yo 

% dw/dteta=non1 dc/dteta=non2 dw/dr=non3 dc/dr=non4 (in phy) 
Yo 

% transpose matrices for column operations in teta direction 
r=sqit(rsquare)’; 

Wey e=C 

non 1=operat(w); 

non2=operat(c); 

wph=phy(w);cph=phy(c); 

non3=[M1*wph'}’; 
non4=[M1*cph']';non4(:,M+1)=-0.5.*sin(teta); 


term 1=-1./r.*non3.*(non2+c_teta); 
term2=1./r.*non1.*(non4+c_r); 
term3=-1.*cotteta./r.*wph.*(non4+c_r); 
term4=-1.*cotteta./r.*non3.*(c_+cph); 
term5=wph./rsquare'.*(non2+c_teta); 
term6=(c_+cph)./rsquare’.*non 1; 
term0=term1+term2+term3+term4+term5+term6; 


ter=[term0;zeros(1,M+1);-1.*flipud(term0(2:N,:))]; 
RR2=fft(ter); 
R2=imag(RR2(1:N,:));R2=R2’; 
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Jo To To To To Yo Yo Yo Yo 


function y=operat(x) 

global N M kmat 

re=zeros(N,M+1); 

X=reti.*x; 

Y=1.*kmat.*X; 
Y_ful=[Y;zeros(1,M+1);flipud(¥ (2:N,:))]; 
y=ifft(Y_ful); 

y=real(y(1:N,:)); 


ToTo To ToT To % To Vo 

function y=phy(x) 

global NM 

re=zeros(size(x)); 

X=re+1.*x; 
Y_ful=[X;zeros(1,M+1);conj(flipud(X(2:N,:)))]; 
dum=ifft(Y_ful); 

y=real(dum(1:N,;:)); 


Jo To To To To To Yo 


function y=solsvd1(u,s,v,b) 
w=diag(s); 
wmin=max(w)* le-12; 
for 1=1:length(w) 
if wi)<wmin;ww(1)=0;else;ww()=1/w();end 
end 
www=diag(ww); 
y=v*www*(u'*b); 


To To To To To Yo Vo 


function y=solsvd(u,s,v,b) 

w=diag(s); 

JYwmin=max(w)* le-8; 

for 1=1] :length(w) 

% if wi)<wmin;ww(i)=0;else;ww(i)=1/w(i);end 
ww(i)=1/w(i); 

end 

www=diag(ww); 

y=v*www*(u'*b); 
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